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Modeling the velocity gradient tensor A — Vu along Lagrangian trajectories in turbulent flow 
requires closures for the pressure Hessian and viscous Laplacian of A. Based on an Eulerian- 
Lagrangian change of variables and the so-called Recent Fluid Deformation closure, such models 
were proposed recently (Chevillard & Meneveau, Phys. Rev. Lett. 97, 174501 (2006)). The result- 
ing stochastic model was shown to reproduce many geometric and anomalous scaling properties of 
turbulence. In this work, direct comparisons between model predictions and Direct Numerical Sim- 
ulation (DNS) data are presented. First, statistical properties of A are described using conditional 
averages of strain skewness, enstrophy production, energy transfer and vorticity alignments, condi- 
tioned upon invariants of the velocity gradient. These conditionally averaged quantities are found 
to be described accurately by the stochastic model. More detailed comparisons that focus directly 
on the terms being modeled in the closures are also presented. Specifically, conditional statistics 
associated with the pressure Hessian and the viscous Laplacian are measured from the model and 
are compared with DNS. Good agreement is found in strain-dominated regions. However, some 
features of the pressure Hessian linked to rotation dominated regions are not reproduced accurately 
by the model. Geometric properties such as vorticity alignment with respect to principal axes of the 
pressure Hessian are mostly predicted well. In particular, the model predicts that an eigenvector of 
the rate-of-strain will be also an eigenvector of the pressure Hessian, in accord with basic proper- 
ties of the Euler equations. The analysis identifies under what conditions the Eulerian-Lagrangian 
change of variables with the Recent Fluid Deformation closure works well, and in which fiow regimes 
it requires further improvements. 

PACS numbers: 02.50.Fz, 47.53. -)-n, 47.27.Gs 



I. INTRODUCTION 

Fundamental understanding of universal features of 
the small-scale structure of turbulence has been a long- 
standing challenge in turbulence research [l|, H, d, 0, H, M, 
While considerable phenomenological understanding 
has been accumulated in recent decades, the challenge of 
relating observed phenomena and statistical properties to 
the dynamical equations (Navier-Stokes) remains unmet. 
The velocity gradient tensor Aij ~ dui / dxj (where u de- 
notes the velocity vector) provides a rich characterization 
of the topological and statistical properties of the fine- 
scale structures in turbulence. Having a spectral peak 
at around the Kolmogorov wavelength fc^ ^ 77""^ (77 is 
the Kolmogorov dissipative length scale), Aij is a quan- 
tity dominated by motions in the viscous range. The 
antisymmetric part of the tensor is the vorticity repre- 
senting small-scale rotation of fluid elements, while its 
symmetric part, the strain-rate tensor, represents fluid 
deformation rate. The Lagrangian evolution of this ten- 
sor can be described by an evolution equation that is 
obtained from taking the gradient of the Navier-Stokes 
equations. The resulting system is unclosed since it con- 
tains the anisotropic part of the pressure Hessian and the 
viscous term. When both of these are neglected, the sys- 



tem is closed (called Restricted Euler - RE- dynamics) 
[1, M, [13, Hii • The RE equations already predict several 
kno wn g eometric turbulence phenomena associated with 
A^j [11 [H, M, M, M, [13, [ll, H [13: such as prefer- 
ential alignments of vorticity with strain-rate eigenvec- 
tors and preponderance of axisymmetric expansion and 
positiveness of the intermediate eigenvalue of the strain- 
rate. Neverteless, RE produces singularities in a rela- 
tively short; finite time (see Ref. [3l for a review). 

Phenomenologically, other phenomena such as small- 
scale intcrmittency may also be probed, by studying the 
probability distribution functions (PDF) of individual 
velocity gradient elements. For instance, it is known 
that the PDFs of longitudinal and transverse gradients, 
e.g. All and A12 in particular directions, respectively, 
can be described by elongated stretched exponential tails 
[1, 121I. [2^ . [23| or by superposition of stretched exponen- 
tial |24l]. Also the moments of these gradients scale in 
non-trivial (anomalous) ways with the Reynolds number 
[25! . [26j . Understanding such anomalous scaling behavior 
of turbulence is another open challenge. Therefore, the 
wealth of geometric, dynamical and statistical turbulence 
phenomena that can be described by the velocity gradient 
tensor, coupled with the fact that a dynamical equation 
(even though unclosed) is available from the gradient of 
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the Navicr-Stokes equations, makes Aij a tensor variable 
of considerable interest for further study. The role of 
pressure in the intermittent nature of velocity gradients 
was also pointed out by Kraichnan in early works [13, HI] . 

Based on prior works [13, M, [H, HO, [Mj, [H, [H, a 
stochastic dynamical model for the time evolution of Aij 
has been proposed [H, The model includes a clo- 
sure for the pressure Hessian, i.e. d^p/dxidxj, and the 
viscous term, i.e. vV^A in terms of the local value of 
the velocity gradient tensor. The approach, reviewed 
in detail in consists in a change of variables from 
Eulerian positions to Lagrangian labels before assuming 
isotropy in the associated gradient tensors to be mod- 
eled. The Eulerian-Lagrangian transformation involves a 
Jacobian matrix that is modeled using the local value of 
the velocity gradient tensor by using the "Recent Fluid 
Deformation" closure. The model system is forced using 
a Gaussian white-in-time random force. The resulting 
stochastic model consists of eight independent coupled 
stochastic differential equations (SDEs) that aim to de- 
scribe the time evolution of each of the tensor elements 
of Aij^ following a fluid particle in a turbulent flow. 

The results of Ref. [U show that the finite-time di- 
vergence exhibited by the Restricted Eulcr system is reg- 
ularized with the inclusion of the proposed models for 
pressure Hessian and viscous term. And, with the ran- 
dom forcing, stationary statistics of the velocity gradi- 
ent tensor are obtained, with realistic statistical proper- 
ties such as preferential alignment of vorticity and the 
preferential state of axisymmetric expansion. The shape 
of probability distribution function of longitudinal and 
transverse gradients is quite realistic and even some well- 
known properties of anomalous scaling in turbulence are 
reproduced [13, [IHl . A limitation of the model is that at 
high Reynolds numbers the resulting distribution func- 
tions became increasingly unrealistic. One approach to 
remedy this problem has been explored [36| by construct- 
ing a multi-scale version of the model, i.e. a matrix shell 
model that describes the velocity gradient tensors in var- 
ious shells at different scales. The closure for the inter- 
scale interaction terms is based on the criterion that the 
total kinetic energy must be preserved by the modeled 
quadratic inter-scale interaction terms. The introduction 
of non-local (in scale) interactions leads to a structure 
of the model that is more difficult to analyze theoreti- 
cally, but it provides an interesting connection between 
the gradients' evolution at various scales and the energy 
cascade mechanism. While the matrix shell model suc- 
cessfully eliminates the problems at high Reynolds num- 
bers, it does not make an explicit connection with the 
physics of the pressure Hessian. At this stage, then, it is 
of interest to further improve our understanding of the 
fundamental properties of the closures proposed in [sj, 
developed from the expression of pressure Hessian and 
viscous term as given from the Navier-Stokes equations, 
within the range of Reynolds numbers in which the model 
of [U works well. 

In i jlHI various model predictions of statistical and ge- 



ometric properties of A beyond those already studied 
in [3^ are compared with Direct Numerical Simulation 
(DNS) at a moderate Reynolds number {TZx ~ 150). 
The model is evaluated using statistical measures already 
studied in [s^l in the context of the "tetrad model". 
These measures include conditional averages of the "dis- 
sipation" (|Sp — SijSij, where S is the symmetric part 
of A) and of the "enstrophy" (i.e. |r2p = ^1^17^^, where 
rj is the antisymmetric part of A). Then, a similar anal- 
ysis is performed with the enstrophy production and the 
strain skewness. The conditional averages are expressed 
in terms of the two principal invariants of A, namely 
R = -(l/3)Tr[A^] and Q = -{1/2)Tt[A^]. Different 
regions in the "(i?, Q)-plane" have distinct physical in- 
terpretations [llj, [ij, 0, [13] and the behavior of the 
computed conditional averages in these different regions 
thus provide useful and statistically meaningful insights 
into the performance of the model in dynamically very 
different regions of the flow. 

In order to quantify and understand the average local 
evolution of the turbulence dynamics in the {R, Q)-plane, 
the probability current of the joint probability density 
V{QtR) is also studied in ijlVl These statistics depend 
explicitly on both pressure Hessian and viscous Lapla- 
cian, and thus the effect of the proposed closures for 
these terms may be compared with the real effects ob- 
tained from the DNS. 

In ^jV] the preferential alignement of vorticity with 
eigendirections of both pressure Hessian and the sym- 
metric part of the viscous Laplacian are studied in de- 
tail. Connections with theoretical results pertaining to 
the Euler equations are also made. Finally, in WII the 
results are summarized and conclusions are presented. 



II. THEORETICAL BACKGROUND 

A. Lagrangian description of the velocity gradient 
tensor 

A description of small-scale structure of turbulence 
based on the velocity gradient tensor Aij = dui/dxj be- 
gins by taking the gradient of the Navier-Stokes equation. 
One then obtains the system: 



dA^ 
dt 



= -A,kAkj - 



dxi dx 



92 Ay 
/ — 

dxkdxk 



(1) 



where d/dt stands for the Lagrangian material derivative 
(i.e. d/dt = d/dt + Ukd/dxk), p the pressure divided 
by the density of the fluid and the kinematic viscosity. 
Because of incompressibility, A must remain trace-free, 
i.e. An = 0. Equation [1] is not closed in terms of A at 
the position x and time t. This can be easily seen not- 
ing that the pressure field is the solution of the Poisson 
equation TT[d'^p/dxidxj] = V^p = —AikAki which shows 
that pressure is highly non-local. And the viscous term 
requires the Laplacian of A which is not known simply 
in terms of A. 



3 



As already mentioned in ^ neglecting pressure Hes- 
sian anisotropy and viscous effects leads to finite time 
singularities because of the strong and unopposed effects 
of the self-streching term — A^. One can find in the 
literature several efforts at regularizing this finite time 
divergence, while keeping the exact self-streching term. 
Firstly, Girimaji and Pope [2^ succcded to do so by con- 
structing a stochastic model with an imposed constraint. 
This constraint is imposed by modifying the non-linear 
term so that the pseudo-dissipation ip = AijA^j [3^ is 
lognormal with a prescribed mean and variance. Inter- 
mittency trends are put in explicitly, by prescribing a 
known variance of log((/j) as function of the Reynolds 
number. 

More recently, two groups proposed the idea that the 
local geometry of the accumulated fluid deformation, 
i.e. formally the Cauchy-Grecn tensor, may represent 
the missing information which allows to regularize the 
RE divergence. Accumulated fluid deformation thus 
has been used to model the pressure Hessian in the so- 
called "tetrad model" of Chertkov, Pumir and Shraiman 
3C, [1^ H^. A similar idea from Jeong and Girimaji 
3l| has been used to model the viscous part of Eq. H]), 
explicitly using the Cauchy-Grecn tensor. Whether or 
not the finite time divergence is regularized, the direct 
use of the Cauchy-Green tensor is limited by the fact 
that this tensor is fundamentaly non-stationary, i.e. as 
time evolves it maintains memory of the initial condi- 
tion. Hence the resultant models for pressure Hessian 
and viscous term are intrinsically non-stationary and de- 
pend on the initial condition chosen to initialize the ma- 
terial deformation tracking. In the following sections, we 
discuss these issues in more detail and also review the 
simplified model of [U that avoids these problems of 
non-stationarity. 

We also point out an alternative approach [4l[ that 
renormalizes the time variable in RE dynamics so as to 
relegate the finite time singularities to infinite time. 



B. Lagrangian mapping and Cauchy-Green Tensor 

Following Refs. [1, [13, S^, one may define a mapping 
Ttg^t between Eulcrian and Lagrangian coordinates: 



X e 



(2) 



where x(X,f) denotes the position at a time t of a fluid 
particle which was at the position x(X, to) = X at the 
initial time to- Given the initial position of a fluid par- 
ticle, this mapping (Eq. is fully defined at any time 

by 



= u(x,t) 



(3) 



A quantity of much interest in continuous mechanics is 
the deformation gradient tensor D, defined as Dij = 
dxi/dXj, which relates the variation of the position of 



a particle when one slightly changes the initial position. 
Differentiating Eq. ^ with respect to Xj, one gets the 
time evolution of D, i.e. 



(4) 



and one can show (3,111] that the Jacobian of the mapping 
7^(,_t, i.e. det(D(t)), is equal to unity at any time by 
virtue of incompressibility, stating that this mapping is 
always invertible. Eq. (|4]) can be exactly solved using the 
product integral |44l| or alternatively, the time-ordered 
exponential |45l. l46l| 



D = [| e''"^(") = T+ exp / dsA{s 



(5) 



The Cauchy-Green tensor, C{t), is defined as the sym- 
metric tensor C = DD^ and its eigenvalue and eigenvec- 
tor system describes the rotation and deformation of ini- 
tially isotropic-shaped fluid volumes into various shapes 
as time goes on. The transport equation of the Cauchy- 
Green tensor can be obtained in straightforward fashion 
113 from Eq. Q: 



C(t)A^(0 



(6) 



Based on the properties of C, studies of isotropic and 
homogeneous turbulence in both numerical [48j and lab- 
oratory 0, [3 flows have shown that cigar (one large 
and two small eigenvalues of C) and pancake (two large 
and one small eigenvalue) shapes are the most common 
shapes of fluid deformation. 



C. Fluid deformation, and pressure Hessian models 

Let us first remark that the pressure Hessian is not 
among the most studied objects in the turbulence litera- 
ture (although, see Ref. [11]). One reason perhaps is that 
it cannot be described naturally from a standard trans- 
port equation along a Lagrangian trajectory. Instead, 
the pressure Hessian is related to the spatial distribution 
of the velocity gradient using singular integral operators 
[H, 50, 51J: 



dxidxj 



-Tr(A2)^ - RV. J %(x - y)Tr(A2)(y)dy 

(7) 

where the integral is understood as a Cauchy principal 
value (P.V.) and fcy is the Hessian of the Laplacian's 
Greens function, namely 



kij (x) 



1 



dxidxj 4tt\x\ 



47r x 



(8) 



One can see from Eq. ([7]) that only the isotropic part 
of the pressure Hessian is purely local (the first term of 
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the RHS of Eq. ([7])). All the nonlocal effects of pressure 
Hessian enter through the anisotropic part (or deviatoric 
part corresponding to the second term of the RHS of Eq. 
([7])). Hence, in this view, the RE approximation can be 
understood as the neglect of all the nonlocal effects im- 
plied by the incompressiblity condition: the correspond- 
ing Lagrangian particle evolves with the flow completely 
independent from its neighbors. As far as we know, the 
tetrad model [l^l is the first model to have been proposed 
for the anisotropic (i.e. nonlocal) part of the pressure 
Hessian. While the authors introduced the model using 
the language of multipoint dispersion of particles that de- 
fine an evolving tetrad shape, a simple interpretation of 
the model can also be given in terms of the deformation 
and Cauchy-Green tensors. 

To begin, one can re-express various Eulerian quanti- 
ties such as the pressure Hessian and the viscous term 
in terms of Lagrangian coordinates, i.e. in terms of the 
fluid particle's position at some initial time to, X. For 
the Hessian tensor of the pressure at the current point 
and time (x, t) one may write 



9^p(x,t) dXp dXg d^pjx, t) 
dxidxj dxi dxj dXpdXq 



d^Xq dpj^.t) 



dxida 



dX„ 



(9) 



The second term entering in the RHS of Eq. ([9|) re- 
quires the knowledge of the spatial distribution of the (in- 
verse) deformation gradient, through its spatial deriva- 
tive. As will be seen later, the adopted approach ne- 
glects short-time variations in the velocity gradient and 
in the context of the proposed Lagrangian model, it is 
natural to neglect spatial fluctuations of the deformation 
gradient, i.e. d^Xq/dxidxj ~ 0. Next, we discuss the 
remaining term of the RHS of Eq. ^ . The fourth-order 
tensor diXpdjXq can be solved along its trajectory using 
the dynamical evolution for the deformation tensors (Eq. 
(|^). For the remaining factor, the Lagrangian pressure 
Hessian d^p/dXpdXg, we choose the simplest assump- 
tion, namely the isotropic assumption: 



d'^p 



1 



d^p 



dXpdXq SdXradX^ 



(10) 



Physically, this assumption states that as time pro- 
gresses, one looses memory about the relative orienta- 
tions of the initial locations X as far as the present value 
of pressure is concerned. The contraction between Spq 
and diXpdjXq then connects the model to the Cauchy- 
Green tensor introduced in the preceding section. 

So far the pressure Hessian can then be rewritten, using 
Eq. (fTO|) , according to 



d^p dXra dX„ c)2p(x, t) 



dxidxj dxi dxj dXmdXn 



= a 



_il d^p{ii, t) 



-idXkdXk^^^ 

wc follow Ref. (3o| 
and use the Poisson equation V^p = —Anm-^mn = 
{l/3)C-^d^p/dXkdXk. Replacing back into Eq. 



To determine d'^p/dX^dX^ 



leads to HO, HI: 

d^p{t) Tr(A- 



dxidxi 



Tr(C" 



2Q 
Tr(C-i) 



(12) 



D. Fluid deformation and modeling the viscous 
Laplacian 

In a similar fashion, following Ref. [sH, this procedure 
can be applied to the viscous Laplacian entering in the 
gradient of the Navier-Stokes equation (Eq. i.e. 



d^A dXp dX„ 

I 7ii — — 

dxkdxk dxk dxk 



52 A 
' dXpdXq 



(13) 



The resulting Lagrangian Hessian of A entering in 
Eq. p3)l will be considered as (i) isotropic, i.e. 
d^A/{dXpdXq) = d^A/{dX^dX„-,)5pj3 and (ii), 
its trace will be modeled by a friction term. i.e. 
d'^A/{dXmdXm) = — 1/^^A. The characteristic length 
scale £ reflects the typical length in the Lagrangian frame 
over which A is correlated. To estimate this length-scale, 
we note that the typical decorrelation time of A along 
its Lagrangian trajectory is known to be on the order of 
tr = (y I^Y^'^ 1 the Kolmogorov time-scale (where e is the 
dissipation rate) 0. During that time, a fluid particle is 
advected by the turbulence over a distance of the order of 
I = u'tk — A, where u' is the root mean square velocity 
(chosen as advective velocity scale) and A the Taylor mi- 
croscale. Finally, recognizing that vj}? = T~^, where T 
is the integral time scale, one then obtains the following 
model for the viscous term: 



d^A 

I 

dxkdxk 



1 Tr(C-^) 



T 



(14) 



This model is similar to the one obtained by Jeong and 
Girimaji [3l| but using a different, more physically mo- 
tivated time scale. 



E. Stochastic model based on the Recent Fluid 
Deformation (RFD) closure 

The various terms entering in the rhs. of Eqs. ([T^ and 
p4)) include the tensor C. If this tensor is obtained from 
its transport equation (Eq. ([5])) subject to the natural 
initial condition Cy (<o) = <^ij, then the closures for pres- 
sure Hessian and viscous term depend strongly on the 
initial time tg, or equivalently. on the initial position X. 
Due to the dispersive nature of turbulent flow, C contin- 
ues to evolve with exponentially growing and decreasing 
eigenvalues. Instead of solving for C from its transport 
equation and having to deal with the problems associ- 
ated with non-stationarity, in (3^ a simple closure was 
proposed. It consists of a sort of 'Markovianization' of 
the dynamics of C in that it is assumed that C evolves 
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in a frozen velocity gradient tensor field during a char- 
acteristic (short) time r. The value of A during that 
time is taken as the most recent value (i.e. the current, 
local, value). And the time-scale chosen is the typical de- 
correlation time-scale of A during its Lagrangian evolu- 
tion, which is known to be of the order of the Kolmogorov 
time-scale, tk- Thus, the initial time is taken to be at 
to = t — Tk, which allows to write in a simple way the 
time ordered exponential entering in Eq. ([5]). We thus 
replace the true Cauchy-Green tensor by a new tensor, 
called the "recent Cauchy-Green tensor' that can be 
expressed in terms of simple matrix exponentials: 



(15) 



This leads to an explicit A-depcndent model for the full 
pressure Hessian: 



dxidxj 



Tr(A^) 
Tr(C;^i) 



and for the viscous Laplacian 

1 Tr(C-i) 



dxi-dxh 



T 



(16) 



(17) 



Inserting Eqs. p6)) and p4|) into Eq. ([T]) and writing 
the equation in the Ito's language of stochastic differen- 
tial equations [H^l the full model for the time evolution 
of the velocity gradient reads 



dA 



Tr(A") 
Tr(CrJ) 



Tr(C 



In 

Tk I 



(18) 

The stochastic time evolution of the velocity gradient ten- 
sor A (Eq. [T51) . as proposed in Refs. [SJ, [11], relates 
the joint deterministic action of the self-streching term 
— A^, the pressure Hessian (Eq. [T6|) and the viscous term 
(Eq. [T7)) . Moreover, the system is forced with a stochas- 
tic Gaussian noise. The deterministic part provides two 
time scales: a small time scale tk and a large one T . 
The latter arises in modeling the viscous diffusion term 
when combining the viscosity with the Taylor-microscalc, 
which in turn is related to the large-scale velocity rms. 
Hence, the deterministic part gives the dependence on 
the Reynolds number 7?,e to the model through the ratio 
{T/tkY ~ T^e, according to classical Kolmogorov dimen- 
sional arguments Q . Dependence on the Reynolds num- 
ber of higher order moments of velocity derivatives (i.e. 
anomalous scalings and the intcrmittency phenomenon) 
have been studied and quantified in Rcf. (35j . The pur- 
pose of this article is to focus on a single Reynolds num- 
ber and to compare it with a DNS fiow (see next para- 
graph). 

The term W is a tensorial delta-correlated noise term 
that has been added in order to represent possible forc- 
ing effects, e.g. from neighboring eddies [sj, UHl- In Ap- 
pendix [XI we describe this noise extensively, and propose 
a way to simulate it. 



F. DNS data and comparisons with the model 

In the following, we will make extensive use of a stan- 
dard direct numerical simulation (DNS) of the Navier- 
Stokes equation for a Taylor based Reynolds number of 
order TZx = 150. Pseudo-spectral simulations are per- 
formed, of an isotropic turbulent flow in a [0, 2ttY box 
using 256"^ nodes. Fourier modes in shells with |k| < 2 are 
forced by a term added to the Navier-Stokes equations, 
which provides constant energy injection rate £/ = 0.1. 
The viscosity of the fluid \s v = 0.00113. The time step 
Ai is chosen adaptively to ensure the Courant number 
A^itniax/Ax ^ 0.15, where Umax IS the maximum veloc- 
ity and Ax is the grid size. In order to make comparisons 
between DNS data and the model, one has to specify a 
value for the parameter of the model tk- At TZx = 150, 
it has been estimated by Yeung et al. that the ratio 
of the Kolmogorov scale and the integral (i.e. velocity 
correlation time scale) time scale is tk/T k, 0.1. Thus, 
in the following, DNS data will be compared to the model 
run with tk = O.IT. Without loss of generality, the in- 
tegral time scale T will be set to unity. It corresponds 
to set time as units of T. The model as written out 
as in Eq. p8|) is solved numerically, with the parame- 
ter Tk — 0.1, using a second order predictor-corrector 
method (see Ref. (S^l) with a time step of At = 10~^. 
One obtains time-series of all of the components of the 
tensor A that display temporally stationary statistics. In 
this article, we have worked with a time-series of length 
~ 10^ in units of the integral time scale T. These can 
then be directly compared to DNS results. Furthermore 
the model provides statistically stationary time-series for 
both pressure Hessian and viscous Laplacian that can be 
also directly compared to DNS data. 



III. CONDITIONAL STATISTICS OF THE 
VELOCITY GRADIENT TENSOR 

We are interested here in studying various properties of 
the velocity gradient tensor conditioned upon the value 
of its two invariants R and Q defined earlier. The joint 
probability density of (R,Q) (the RQ -p\a,ne) has been 
studied in the past [l^, [ll|, Tisl . [Tsl . [STj and can be used 
to characterize the frequency of occurrence of the various 
local topologies of the flow. 

For instance, in a simple way, the second invariant 

= -iTr(A2) = i|c.p-iTr(S2) (19) 

can be understood as the competition between enstrophy 
{uj denotes vorticity) and dissipation (per unit viscosity). 
Then, positive Q represents rotation-dominated regions 
and negative Q dissipation-dominated regions. Analo- 
gously, the third invariant 



R 



1 



Tr(A3) 



1 



1 



Tr(S-'') 



(20) 
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FIG. 1: Joint PDF P(Q*,i?*) of R* = R/{SijSij}^^^ and 
Q* ~ Q/{SijSij) calculated from DNS (a) and the present 
model (b). Contour lines are the same in the two cases, log- 
arithmically spaced by a factor of 10, and start at 10 near 
the origin. The thick line represents the zero-discriminant (or 



Vieillefosse) line: ^R^ + Q-^ 



■ 0. 



represents competition between enstrophy production, en- 
tering in the enstrophy evolution (tI. [s^. i.e. 



ld\iu\'- 
2 dt 



— OJiSijLOj 



(21) 



and the dissipation production (or Strain skcwness |30l |) 
entering in the dissipation evolution , i.e. 



dTr(S2 
dT 



3 1 

-2Tr(S )~-ujiSijUjj^ 



'2Si' 



d'^p 
dxidxi 



-vSijV Sij 



(22) 

Let us remark that one may interpret the RQ plane in a 
different way, based on the eigenvalues of A (two of them 
can be the complex conjugate) and the zero-discriminant 
line (i.e. the "Vieillefosse" line, namely ^i?^ + = 0). 
Sec for instance [ll|, [H, [13] . 



A. The joint PDF in the RQ-plane 

We show in Fig. [1] the joint PDF of R and Q, 
or equivalcntly the joint PDF 'P{Q*,R*) of the non- 
dimensionalized invariants R* ~ R/ {SijSij)'^^'^ and Q* ~ 
Q/{Si]S^j), for both DNS and the model (Eq. dUl)). In 
order to compute various conditional averages, a range 
[— 1 ; 1] of values of the two non-dimensionalized invari- 
ants R* and Q* of A is discretized in 25 equally spaced 
bins. For the DNS, one observes the predominance of 
the Enstrophy-Enstrophy production quadrant {R* < 
and Q* > 0) and the Dissipation-Dissipation produc- 
tion quadrant {R* > and Q* < 0). The predomi- 
nance of these two quadrants has been observed before 
in the literature [13] ■ The model reproduces these ba- 
sic trends fairly accurately, with the characteristic "tear- 
shape" elongation along the "Vieillefosse tail" in the 
Dissipation-Dissipation production quadrant. But one 
also observes that the model overestimates the total prob- 
ability in the Enstrophy-Dissipation production region 
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FIG. 2; Isocontours of the conditional enstrophy 
{uJiUJi\Q* , R*)V{Q* , R*) where oj is the vorticity, and 
conditional dissipation {SijSij\Q* , R'')'P{Q* , R*). Following 
Ref. [30| . both quantities are renormalized by (2flij^ij), 
where ft is the rate of rotation tensor. Level of contour lines 
are 0.15, 0.3, 0.45, 1, 2, 3, and 4. 



(i.e. R* > and Q* > 0) and underestimates the 
Dissipation-Enstrophy production region (i.e. R* < 
and Q* < 0). It will be shown later that this is caused 
by limitations in how the pressure Hessian is closed and 
modeled. 

Next, we check whether or not dissipation (resp. en- 
strophy) is dominantly associated with the R* > 
and Q* < (resp. R* < and Q* > 0) quad- 
rants. Following the approach already used in (soj 
we present in Fig. [5] the conditional averages of dis- 
sipation, i.e. {Tt{S'^)\Q*,R*)V{Q*,R*) and enstrophy 
{LLiiUJi\Q*,R*)V{Q*,R*) on R* and Q* . Averages are 
weighted by the joint density V{Q*,R*) to ensure that 
the sum over all possible values of R* and Q* gives the 
averages of respectively dissipation and enstrophy. We 
clearly see that the quadrants i?* > 0; Q* < (or R* < 0; 
Q* > 0) are dominated by dissipation (or enstrophy, re- 
spectively). The model reproduces these conditional av- 
erages quite accurately. 



B. Enstrophy production, Strain Skewness and 
Energy transfer 

A similar study is performed with the various quanti- 
ties entering in the definition of the third invariant R (Eq. 
(j20p ). namely the enstrophy production and the strain 
skewness [3^. In Fig. [3] (a) to (d) these various quan- 
tities are shown, together with the predictions from the 
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FIG. 3: Isocontours of the conditional Strain Skew- 
ness -(Tr(S^)IQ*, J?*)-p(Q*,_R*), enstrophy produc- 
tion {ujiSijUjjlQ* , R*)'P{Q* , R*) and energy transfer 
~{Tr{A'^A'^)\Q,R)ViQ*,R*). Following Ref. var- 
ious quantities are normalized by the average transfer 
|(TV(A^A^)>|. 

stochastic model. In all cases, it is apparent that model 
predictions are quite accurate. We see also that enstro- 
phy production is clearly dominant in the R* < and 
Q* > quadrant. Let us mention that in the cnstrophy- 
dissipation production dominated region (i?* > and 
Q* > 0), entrophy production becomes weakly nega- 
tive, stating that in this region, enstrophy decreases with 
time (see Eq. ([?!]) ). Also, in Fig. [n](a-b), we see that 
strain skewness is dominating in the bottom-right quad- 
rant, but remains very important in the top-left quad- 
rant. This is mainly linked to the fact that the evolution 
of dissipation not only depends on the strain skewness 
(or dissipation production), but also on enstrophy pro- 
duction, and a term linked to the pressure Hessian (see 
Eq. (HID). 

A related quantity of interest is 

- Tr(A2A^) = -Tr(S''') - ^uoiS^jOJj (23) 

which describes the time evolution of the pscudodissipa- 
tion dTr(AA^)/dt = -Tr(A^A^) in the RE approxima- 



tion. This quantity is sometimes called "energy transfer" 
jsol . [ssl . [56| when A is defined by filtering in the inertial 
range in the context of large eddy simulations (see [57l|). 
While here A is not filtered and therefore no such direct 
physical interpretation is available, this quantity is still 
presented as additional documentation of the properties 
of A. Results are displayed in Fig. [3]Je-f). Once again, 
the model reproduces well the trends observed in DNS, 
including negative regions in the top- right quadrant. 



C. Geometric alignments of vorticity with 
strain-rate eigenvectors 

An important universal feature of fully developed tur- 
bulent flows is the preferential alignment of vorticity 
along the eigendirection of the intermediate eigenvalue of 
the strain-rate tensor S (see Q and references therein). 
To study the alignment properties of vorticity condi- 
tioned on various values of R and Q the (i?, Q) plane is 
divided into four regions related to the eigenvalue struc- 
ture of A. Instead of the (5 = line to separate high 
and low rotation regions as was done in the qualitative 
discussions of the previous sections, we now use the quan- 
titatively more precise classification, in which the (i?, Q) 
plane is divided into high and low rotation regions by the 
zero-discriminant line, i.e. Q = — (27i?^/4)^/'^. 

Fig. |4] shows the PDF of the cosine of the angle be- 
tween vorticity and eigendirections with the most neg- 
ative (a-b), intermediate (c-d) and most positive (e-f) 
eigenvalue of the stress, for both DNS and the model. 
The different symbols denote the results obtained in sep- 
arate quadrants as separated by the Vieillefosse (i.e. zero 
discriminant) and the R = lines. In Fig. IHJg-h) is dis- 
played the unconditional PDF independent on the quad- 
rant, i.e. as obtained in all regions. As already ob- 
served in [sl], the model predicts accurately the pref- 
erential alignment with the intermediate eigendirection 
(solid line), a trend of being orthogonal to the most con- 
tracting direction (dashed line), and an almost entirely 
dccorrelated trend with the most extensive eigendirection 
(dash-dotted line). The agreement between DNS and the 
model is excellent in all cases, even when conditioning on 
the separate quadrants. It is interesting to note that in 
(a) and (b), as well as in (e) and (f), the alignment PDF is 
essentially the same in three quadrants but very different 
in the i? > and Q > -(27i?V4)^/^ quadrant. In Fig. 
mja) and (b) we observe that while the vorticity is mostly 
perpendicular to the most contracting eigendirection, in 
the top-right quadrant the vorticity is in fact not orthog- 
onal to the contracting cigen-direction. This is the "vor- 
tex contracting" quadrant with an unstable focus and 
one contracting direction. This would suggest that the 
vorticity is aligned with the contracting direction. There 
is instead no strong preferred alignment but there is an 
almost zero probability that the vorticity is perpendic- 
ular to the contracting eigendirection. But, on average, 
when taking into account all the possible values for R and 
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FIG. 4: PDF of the cosine of the angle 6 between vorticity 
and the different eigendirections of S: (a) and (b) for the most 
contractive eigendirection (negative eigenvalue), (c) and (d) 
for the intermediate eigendirection, and (e) and (f) for the 
most positive (extensive) eigen-direction. As it is schemat- 
ically displayed in the inset of Fig. Hfb), different symbols 
are obtained from the four different regions of the {R, Q) 
plane delimited by the i? > and the Vieillefosse (or zero- 



discriminant) lines given by Q 



27^2^,1/3, 



{R> and 



Q > -(f i?^)'/^), V (i? > and Q < -(f 7?^)'/=*), □ (J? < 
and Q > -(f i?2-,i/3)^ (R<0 and Q < -(f i?2)i/3) 
and (h) show the unconditional PDF over the entire {R, Q) 
plane. Different lines correspond to different associated eigen- 
values: most negative (dashed), intermediate (solid) and most 
positive (dashed-dotted). 



Q; (Fig- lllg-h))i vorticity remains weakly orthogonal to 



this eigendirection. In terms of the alignments with the 
intermediate eigendirection, one may have expected the 
preferential alignments to come mainly from the bottom- 
right quadrant as predicted by the asymptotic diverg- 
ing state of the RE equations (lol. [Tl|. Nevertheless, in 
Fig. ^iid (d) we observe instead that the alignment 
with the intermediate eigendirection occurs quite inde- 
pendently of the characteristic values in the (i?, Q) plane. 
In Fig. |4je) and (f) we observe that the alignment with 
the most extensive strain-rate eigendirection is random 
in all quadrant except, again, in the top-right quadrant 
where i? > and Q > -{21R^ / A^f^ . Vortex "contrac- 
tion", when it happens, appears to occur because it is 
mostly orthogonal to the extensive direction and also 'not 
orthogonal' to the contracting direction, rather than be- 
ing preferentially aligned with the contracting direction. 
The stochastic model predicts these non-trivial statistical 
geometric behaviors quite well. 



IV. PRESSURE HESSIAN AND VISCOUS 
TERM 

Let us now focus directly on the terms requiring clo- 
sure, namely the pressure Hessian and the viscous term, 
instead of the statistics of the velocity gradient tensor 
considered in the previous section. One option could be 
to compare individual realizations of the model terms 
with the corresponding DNS values along Lagrangian tra- 
jectories. However, since these terms fluctuate greatly in 
the DNS, a statistically more robust comparison is per- 
formed using conditional averages, conditioned on R and 

g. 



A. Probability current and conditional averages 

The approach used in Ref. [ll] is followed, based on 
a Fokker-Planck equation for the dynamics of R and Q. 
To summarize the approach, we notice that along a La- 
grangian trajectory, appropriately contracting ffl with 
A and , using the Cayley-Hamilton theorem , one 
can show that the time evolution of the invariants R* 
and Q* are given by 



dQ* 
dt* 



-3R* 



1 



1 



and (24) 



dt* 3 ' 



a a 



(25) 



where = (SijSij) is the strain variance and t* = at 
the non-dimensional time. Also, stands for (minus) 
the deviatoric part of the pressure Hessian, i.e. 



\dxidxj 



3 dxhdxi 



(26) 



and H"^ = vV^A. is the viscous term (recall that in the 
RE approximation, = = 0). The Fokker-Planck 
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It can be decomposed into W = Wre + Wp + Wi^, with 



Q\R*)V{Q\R*) (29) 



which describes the deterministic (closed) part of the evo- 
lution of the two invariants, 



Q\R*)V{Q\R*) , (30) 



describing the pressure Hessian effects on the evolution 
of R* and Q* and finally, 



Q*,R* )ViQ*,R*) , (31) 



describing the effects of the viscous term. An additional 
current might be considered in this description, linked to 
an additional forcing term that has been neglected in the 
Navier-Stokes equations (Eq. ([T])). This forcing is indeed 
negligible in front of the other terms of the rhs. of Eq. 
^ since it can be written as the (small-scale) gradient of 
the large-scale forcing of the velocity, and thus, we will 
neglect its associated probability current. 

Conversely, in the Fokker-Planck equation (Eq. ^7} ) 
for the joint probability distribution of R* and Q* ob- 
tained from the model (Eq. (HH])), one has to take 
into account another term which comes from the delta- 
correlated Gaussian forcing. Appendix |B] provides the 
required background needed to compute the probabil- 
ity flux resulting from the stochastic forcing term in our 
model, i.e. the diffusion terms entering in Eqs. (j24p 
and . It is shown that the currents associated to the 
deterministic and random parts of the joint stochastic 
evolution of R and Q predicted by the model (see Eq. 
(|B4p ) are of the same order of magnitude. 



FIG. 5: Vector and streamline plots of the probability current 
associated to the (a-b) RE approximation (Eq. (|29p ). (c-d) 
the pressure Hessian (Eq. (|30|l '). (e-f) the viscous term (Eq. 
^T^). The total current (Eq. (^5)) ') is represented in (g-h). 
The scale of the vectors is the same in all plots and a reference 
is given below Fig. (h), whose (non-dimensional) magnitude 
is 2.10"\ 



equation describing the time evolution of the joint density 
V{Q*,R*) may be written as 

.W - , (27) 

-■OR' ' 

where the divergence of the probability current W con- 
trols time variations of the joint probability density V. 
The probability current can be written in terms of con- 
ditional averages as 



dV 

dt* 



■ \ dt' 



Q\R* )V{Q\R*). 



(28) 



B. Results 

In Fig. O the vector plots and associated streamlines 
corresponding to the various probability flux terms are 
presented. Both results obtained from DNS and from 
the model are shown. 

First, as reference we present in Fig. [5]Ja-b) the closed 
RE current VVre (Eq. ^). As is well-known 
the deterministic Wre probability current pushes prob- 
abilities towards the right tail of the Vieillefosse line. 
Since the model predicts acurately the joint probabil- 
ity V{Q* ,R*), agreement between DNS and the model 
predictions (length of vectors) is quite good because the 
self-streching term — is taken into account exactly in 
the model (Eq. [TH]) . 

The action of the pressure Hessian, given by the prob- 
ability current Wp and shown in Fig. [5jc-d), is quite 
interesting. From the DNS data, two main pressure Hes- 
sian effects can be observed. First, the pressure Hessian 
counteracts the effects induced by the RE terms since the 
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flux goes towards the center of the RQ plane along the 
right tail Vieillefosse line. This feature is well reproduced 
by the model, with vector magnitudes of the same order. 
Another important effect of the pressure Hessian is that 
in the R < left half-plane, the probability current leads 
the probability towards the left tail of the Vieillefosse 
line, namely towards dissipation-cnstrophy production 
dominated region (in lower- left direction). This feature 
is not reproduced by the model, which instead appears to 
act exclusively in vertical direction, upward in the Q < 
plane, and downward in the Q > side. This explains 
perhaps why the model leads to an underestimation (see 
Fig. [T|) of the probability of dissipation-enstrophy pro- 
duction events (i.e. the bottom-left quadrant). At the 
same time, the absence of "left-ward" flux out of the 
enstrophy-dissipation production region (i.e. top-right 
quadrant) may explain why the model ovcrpredicts the 
probability of events in that quadrant. A very marked 
feature of the DNS results is that the magnitudes of the 
vectors are essentially negligible in the entire vortex con- 
traction quadrant above the right Vieillefosse line (i? > 
and Q > — (27i?^/4)^/'^), leading to some uncertainty in 
the computed streamlines there. 

Another main difference between DNS and model pre- 
dictions is the fact that for the model, VVp vanishes at 
vanishing Q, but not for the DNS. A general feature of 
the pressure Hessian model is that its deviatoric part is 
directly proportionnal to Tr(A^) = -~2Q (see Eq. (fT6|) '). 
Incidentally, the same occurs in further generalizations 
that have been proposed by Gibbon and Holm (see 
their Eq. (5.8)), namely 



Wp - 



' N 



Tr(G„) 



Tr(A^) with 



N 
n=l 



(32) 

where the scalars cn are undcrtcrmincd and G„ arc any 
non-singular symmetric tensors. Once again, we can see 
that the deviatoric part of the pressure Hessian is still 
proportional to Tr(A^). For the sake of completeness, let 
us remark that at least formally, this issue does not arise 
in the matrix shell model of [3g . This is because the non- 
local closure terms in the matrix shell model [111 are not 
directly proportional to Tr(A^) since the connection to 
pressure Hessian and Poisson equation for pressure is not 
included in that approach. It would be very interesting 
to check if the comparison with DNS for the equivalent 
probability current in the matrix shell model, i.e. based 
on relevant portions of the quadratic nonlinear interac- 
tion terms, is better or not. Such studies arc left for 
future work. 

To more clearly isolate the behavior of the pressure 
Hessian near Q = line, we study the magnitude of the 
anisotropic (i.e. deviatoric) part of the pressure Hessian, 
conditioned on the local value of Q. Fig. [6] shows the 
conditional average of the norm (square) of the devia- 
toric part of the pressure Hessian, i.e. ^|H^|^ where 

|HP|2 = Tr (hp (HP )^) , as a function of the local value 
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FIG. 6: Conditional average of the norm square of the devia- 
toric part of the pressure Hessian, (|H''|^ \Q\ with respect to 
the invariant Q. Both |H''|^ and Q, are non-dimensionalized 
by their respective standard deviations. DNS results (solid 
line) and model predictions (dashed line) are shown. 



of the invariant Q, for both the DNS and the model. 
For vanishing Q, the conditional average jH^p from the 
DNS docs not vanish. As discussed before this property 
is not reproduced by the existing models for the pressure 
Hessian, namely the tetrad model (Eq. (fT2|) ). the CM06 
model (Eq. (fTB]) ') and the generalized tetrad model (Eq. 
([5^ ) because all of them predict a pressure Hessian pro- 
portional to Q. Secondly, one can see that for the range 
of Q under consideration (i.e. Q € [—eg, erg] (erg stands 
for the standard deviation of Q), the conditional average 
of |H''|^ behaves as ~ \Q\ (up to a positive additive con- 
stant) whereas the model predicts it to be proportional 
to . The fact that the model predicts a quadratic be- 
havior can be understood from a Taylor's development, 

namely « — Qr^^S leading to (^jH^I^ ~ since 

(t^ |S|^) ^ 1. The small asymmetry in the quadratic 
behavior seen in Fig. [H] is caused by higher order terms 
entering in the expansion for the model. 

In terms of the viscous term, one observes in Fig. [5je-f) 
that the model reproduces the probability flux reasonably 
well. Consistent with the observations already made in 
[l6t the viscous effect is to push the probabilities towards 
vanishing R and Q, not only along the Vieillefosse line, 
but everywhere. We notice that the model ovcrpredicts 
the magnitudes, i.e. at this Reynolds number the model 
provides too strong damping but is qualitatively correct. 

In Fig. EKg-h) is shown the sum of all these terms, 
namely the total probability current W ~ re + VVp -I- 
Wiy. For the model case, another term coming from the 
Gaussian delta-correlated forcing (see Appendix |B] and 
Eq. (jB9p ) has been added. The circular motion around 
the origin of the RQ plane has already been reported in 
Refs. [lljllll. At this point, from Fig. [5]one can observe 
that all the terms (self-stretching, pressure Hessian and 
viscous Laplacian) entering in the Navier-Stokcs equa- 
tions (Eq. IH)) arc of the same order of magnitude (vis- 
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cous Laplacian is a little bit smaller than the two other 
terms, but not by much). Similar conclusions can be 
drawn for the deterministic terms entering in the model 
(Eq. (fTS]) and figs. [DJb-d-f)), although, as shown in ap- 
pendix [B] the amplitude of the forcing is not negligible 
either. Focusing on the total probability current [5jg-h), 
we reach the conclusion that the fact that the modeled 
pressure Hessian (fig. [S^d)) is not able to reproduce the 
probability flux towards R* < and Q* < regions as 
it is observed in DNS (fig. EJ^c)), explains why the model 
over predicts R* > and Q* > regions and under pre- 
dicts R* < and Q* < regions, as observed in Fig. [T] 
We also remark that streamlines shown in the left part 
of Fig. \Mjn) are not significant because of the very low 
values of the joint probability V{Q* ,R*). 



V. VORTICITY ALIGNMENTS WITH 
PRESSURE HESSIAN AND VISCOUS 
LAPLACIAN EIGENDIRECTIONS 

A. Pressure Hessian 

Here we focus on vorticity alignment properties along 
the eigendirections of the pressure Hessian. It has been 
derived, in the inviscid limit (Euler equations) [sol. Is^. [60l| 
that vorticity uji tends to be simultaneously an eigenvec- 
tor of the rate of strain tensor S and the pressure Hessian. 
When v = Q (see Eq. ([2T|)). 



DNS 



Model 



dhJi 

~dt 



— SijUJj , 



(33) 



and taking another time derivative and using the time 
evolution of S (sec Eq. (|22|) ). we get 



d LOi 



(34) 



Following Ref. |50| , we then notice that if vorticity of a 
fluid particle continues to be an eigenvector of the rate- 
of-strain tensor, then it is also an eigenvector of the pres- 
sure Hessian. To see if such a trend is observed in a finite 
viscosity turbulent flow we will quantify the aligments of 
vorticity with the eigendirections of the pressure Hes- 
sian. Such an analysis based on DNS has been already 
performed [6l|, [G^I , but here the purpose is to compare 
results with predictions of the model. 

Let us focus on alignment properties of vorticity with 
respect to the eigendirections of the deviatoric part of 
the pressure Hessian (i.e. —Hf^ defined in Eq. (|26p ). 
Alignment PDFs are shown in Fig. [71 presented in a 
similar fashion as in Fig. U) We can see in Fig. [7]^a) 
that vorticity is preferentially orthogonal to the eigendi- 
rection of the smallest eigenvalue, except in the top-left 
quadrant (i.e. i? < and Q > -(27/47?^)^/^) where the 
local topology is dominated by one direction of streching 
and a stable focus. Also in this quadrant, vorticity is 
preferentially aligned with the extending eigendirection. 
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FIG. 7: PDFs of the cosine of the angle between vorticty 
and the different eigendirections of pressure Hessian tensor 
H'': (a) and (b) for the smallest eigenvalue eigendirection, 
(c) and (d) for the intermediate eigendirection, and (e) and 
(f) for the most positive eigen-direction. As it is schemati- 
cally displayed in the inset of Fig. Elb), in a similar fashion 
as in Fig. 2] different symbols are obtained from the four dif- 
ferent regions of the {R, Q) plane delimited by the _R > 
and the Vieillefosse (or zero-discriminant) lines given by 
Q = -(f i?2)i/3. o (i? > and > -{ilR^f/^), V (J? > 
and Q < -(f J?2)i/3), □ (J? < and Q > -(f i?2.,i/3)^ ^ 

(i? < and Q < -(^i?^)^/^). Figs (g) and (h) show the un- 
conditional PDF over the entire {R, Q) plane. Different lines 
correspond to different associated eigenvalues: most nega- 
tive (dashed) , intermediate (solid) and most positive (dashed- 
dotted). 
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The model predicts a slightly different picture since vor- 
ticity is predicted to be also preferentially othogonal to 
the eigendirection except in the top-right quadrant (i.e. 
i? > and Q > -(27/4i?2)i/3) foj- which local topol- 
ogy is dominated by one compressive direction and an 
unstable focus. 

In Fig. IZl^c-d), we focus on the pressure Hessian 
eigendirection of its intermediate eigenvalue. In a sim- 
ilar way as with eigendirections of the rate-of-strain. vor- 
ticity is preferentially aligned with this eigendirection in 
all the quadrants, and this is also very well predicted by 
the model. About the eigendirection of the largest eigen- 
value, we can see that for the DNS (Fig. [7]Je)) vorticity 
is weakly preferentially aligned with the eigendirection, 
except again in the top-left quadrant where vorticity is 
clearly preferentially orthogonal to this eigendirection. In 
Fig. El^f), we can see that the model predicts most of the 
trends quite well in all the quadrants. For the average 
results over all quadrants, it can be seen in Figs. Eljg-h) 
that the model predicts with a fairly good accuracy the 
behavior of vorticity with the eigendirection of the in- 
termediate eigenvalue (although it overpredicts the peak 
a little bit). In the other extremal eigendirections, the 
model reproduces the moderate peak at cos{9) ^ but 
misses the narrow peaks near alignment at cos(0) ~ 1. 

The fact that the model reproduces very well the events 
for which vorticity happens to be an eigenvector of the 
rate-of-strain tensor can be understood phenomenologi- 
cally in the following way. When vorticity is an eigen- 
vector of the rate-of-strain, i.e. So; = /3a;, then vorticity 
is also an eigenvector of the velocity gradient tensor it- 
self, namely Alj = (S + n)Lj = Pu), since by definition 
ftu} = 0. Let us notice that vorticity is also an eigen- 
vector of = S — ri with the same eigenvalue /3. It is 
then straightforward to show by induction that for any 
power n g N, A"u; = (A^)"a; = f3"u}. For these very 
particular events in which vorticity is an eigenvector of 
the rate-of-strain, one notices than the matrix exponen- 
tial entering the (inverse) recent Cauchy-Green tensor (in 
Eq. can be written as 

+°° ( \n+m 

^-i^^-r^A g-TK-A^ L^EL (A^)"A'" 

n\ra\ 

n,m— 

(35) 

From Eq. p5|) it is easily seen that vorticity is also an 
eigenvector of both the recent Cauchy-Green tensor (Eq. 
(fT5|) ) and its inverse. For example, C^^o; — e~^'^'^''^u). 
Finally, since the pressure Hessian is modeled as propor- 
tional to (Eq. HH)), we can state here that the 
present model is such that when the vorticity is an eigen- 
vector of the rate-of-strain, then it is also an eigenvector 
of the modeled pressure Hessian, the ordering of the asso- 
ciated eigenvalues being respected in absolute value. More 
precisely, if vorticity is an eigenvector of S with respec- 
tive eigenvalue then vorticity is also an eigenvector of 

the pressure Hessian with eigenvalue " ^"^"^^"'^^ 
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FIG. 8: PDF of the cosine between vorticity and eigendirec- 
tions of i^AS associated to the smallest (dashed line), inter- 
mediate (solid line) and biggest (dot-dashed line) eigenvalues 
of the strain-rate tensor. 



B. Viscous Term 

Let us now focus on the geometrical properties of the 
viscous tensor, namely ^V^A, appearing in Eq. ([T]). Let 
us begin with its symetric part ^^V^S. The present model 
(i.e. Eq. p8)) ) contains the closure for the viscous term 
pointed of Eq. (fT4|) that stated that the Laplacian of 
A is proportionnal to A itself. In terms of eigendirec- 
tions, it is assumed that both S and V^S have the same 
eigendirections. Among others, alignment properties of 
vorticity lj and eigendirections of V^S should be exactly 
the same as alignments of vorticity with eigendirections 
of S. To determine whether this is observed in DNS flows, 
we present in Fig. [8] PDFs of the cosine of the angle be- 
tween vorticity and eigendirections of the viscous term, 
for both DNS (a) and the model (b). Fig. (Hb) is in 
fact the same as Fig. IH^h) and is reproduced here for 
convenience. We see clearly that the overall geometrical 
picture are really close between DNS and the model, and 
that there is preferential alignment of vorticity with the 
eigendirection associated to the intermediate eigenvalue 
of the Laplacian term as well. 

Let us now focus on the antisymmetric part of the vis- 
cous term, namely j/V^il. The vorticity vector is given 
by uji = —^Eijk^jk- The viscous term -cc = i/W^uj can 
also be written as a vector, i.e. rui = —^Eijk'^^'^^jk- 
The model for the viscous term (Eq. ((T7)) ) implies that 
the angle 9 = (177,0;) between the vorticity and the vor- 
ticity Laplacian is fixed and equals tt, i.e. same direction 
but opposite orientation, since the model (Eq. p7p ) is 
proportional to the velocity gradient tensor A with a 
negative coefficient. One may wonder if this is consis- 
tent with DNS data. We represent in Fig. M the PDF 
of cos 9 estimated from the same DNS fields. Clearly we 
see that the two vectors share preferentially the same di- 
rection, but opposite orientation. The model is therefore 
consistent with the observed alignment trends of the full 
Laplacian of velocity gradient. 
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FIG. 9: PDF of the cosine of the angle 9 between vorticity 
uJi = —^Eijk^jk and tn, the Laplacian of vorticity vector 
vji = —^Eijki'V^fljk, obtained from DNS. 



VI. CONCLUSIONS 

Extensive comparisons have been made between pre- 
dictions of a new stochastic Lagrangian model for the 
velocity gradient tensor and results from DNS at a cor- 
responding moderate Reynolds number. The model re- 
produces many inherent geometric and statistical prop- 
erties of small-scale turbulence quite well. The statistics 
of alignment angles between vorticity and the principal 
axes of the rate-of-strain tensor are very well reproduced. 
The joint statistics of velocity gradient invariants R and 
Q are also reproduced well. Specifically, the joint PDF's 
elongation into the top-left and bottom-right quadrants 
observed in the DNS also occurs in the model. Some dif- 
ferences occur in the model in the top-right and bottom- 
left quadrants. In order to directly assess the action of 
the modeled pressure and viscous terms in a statistically 
robust fashion that takes into account the local topology 
of the flow, the probability current W has been stud- 
ied. The agreement between DNS and model predic- 
tions is good near the dominant Vieillcfosse tail in the 
lower-right quadrant of the {R, Q) plane. However, in 
the dissipation-enstrophy production dominated region 
(bottom- left), the model docs not reproduce the true dy- 
namics and requires further developments. Finally, the 
alignment properties of vorticity with respect to the prin- 
cipal axes of the pressure Hessian tensor have been stud- 
ied. The model reproduces quite well the preferential 
alignment of vorticity with the eigendirection associated 
to the intermediate eigenvalue. We elucidate the fact 
that in the model an eigenvector of the rate-of-strain is 
also an eigenvector of the pressure Hessian, and this is 
in fact consistent with known behavior of vorticity in the 
inviscid limit (i.e. the Euler equations). 

This analysis has confirmed that the stochastic model 
is capable of predicting many non-trivial properties of 
small-scale turbulence as described by the geometric and 
statistical properties of the velocity gradient tensor. Nev- 
ertheless, there appear to be difficulties in specific regions 



of the flow, especially those in which the vorticity is being 
contracted such as in the top-right or bottom-left portion 
of the invariant (i?, Q) plane. Whether these drawbacks 
of the model are also related to the difficulties observed 
when raising the Reynolds number of the fiow also re- 
mains to be explored. 
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APPENDIX A: DEFINITION AND 
IMPLEMENTATION OF TENSOR GAUSSIAN 
FORCING 

In this appendix, the tensorial Gaussian forcing dW 
entering in the model (Eq. p8)) ) is described. It can be 
written as 

dW^j = D.jkidBki (Al) 

where Dij^i are the diffusion coefficients and dB is a ten- 
sorial isotropic Wiener process, whose components are 
such that 

{dB,j) = and {dB.jdBki) = IdtS^kSji ■ (A2) 

The coefficients Dijki arc chosen such that the noise dW 
is consistent with a trace- free, homogeneous and isotropic 
tensor, of a given unit variance, namely {dWijdWki) = 
2dtDijpgDkipq with 

DijpqDkipq — 26ik5ji — -SijSki — -SiiSjk ■ (A3) 

As a conscquency, longitudinal components of the noise 
dW are of variance 2dt, and Adt for the transverse ones. 
Also, let us recall that the dimension of the diffusion co- 
efficients is time~^/^. If the tensor D is assumed isotropic 
itself, then the unique solution of Eq. (|A3[) is given by 

Dzjpq = aSijSpq + bStpSjq + cS^qSjp , (A4) 

with 

^_1 3 + Vl5 VTo + V6^^_ 1 

"~ syio-H V6' ^ 4 ''^^ + ' 

(A5) 



APPENDIX B: STOCHASTIC DIFFERENTIAL 
EQUATIONS 

The basic tensorial stochastic differential equation 
(SDE) Eq. (HH) can be written as 

dA,j ^ Vi-jdt + D.jkidBki , (Bl) 

where Vij is the drift coefficients representing the sclf- 
streching, pressure Hessian and viscous terms entering in 
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FIG. 10: Vector and streamline plots of the probability cur- 
rent for the model associated to the (a) the drift coefficients 
and (b) the diffusion coefficients (see text and Eq. (|B9|l l. 
The scale of the vectors is the same, reference is given in each 
figures, whose (un-dimensionalized) magnitude is 2.10^^. 



To compute the time evolution of the invariants R and 
Q, vsfe use Eq. (|B3p with the particular (time indepen- 
dent) nonlinear transformation — Q = — Tr(A^)/2 = 
-A,,A,,I2 and ii=R^ -Tr(A3)/3 = -Ay A,feAfe,/3. 
We get 



iDijpq 

V- ■ A ■ A ■ 



dt 



AliDijpqD jlpq\ dt 



(B4) 



-^jr Af iDijpq dBpq 



- fjP A A — A- A - 
drift terms coming from 



where, using former notations, we notice that —VijAji = 
—3R— HfjAji — H^jAji and using the Cayley- Hamilton 
theorem, -VijAjqAqi = 
Furthermore, the "spurious" 
the delta-correlated Gaussian noise vanish, i.e. using 
Eq. (|A3p one can show that DijpqDjipq = and 
AiiPii pnDji nn = 0. From a straightforward manner 
[29I . I52 . [SSf . , we get from Eq. ()B4p the corresponding 
Fokker-Planck equation for the joint probability 'P{Q,R) 



dV 
'dt 



d_ 



[PN,] 



d 



(B5) 



In Eq. (jB5|) . the new coefficients Ni and My can be easily 
obtained from Eq. (|B4p . although, they are not known as 
functions of = Q and £,2 = R- Therefore, we will use 
conditional averages to estimate them. Henceforth, we 
will write the Fokker-Planck equation with conditional 
averages and obtain 



Eq. ([T5| and DijkidBki is the forcing term described in 
Appendix [X] The associated Fokker-Planck equation for 
the joint probability of the velocity gradients /(A; t), of 
is given by 

[fVv] + . [fDvpqDkipq] ■ (B2) 



dl 
dt 



a 



OA, 



OA^jdAki 



Here we are interested in the joint probability V{Q,R) 
of two invariants of the velocity gradients, namely R and 
Q (c.f. Eq. (P7)) ). In Ito interpretation, both the stochas- 
tic equations governing the dynamics of R and Q and the 
associated Fokker-Planck equation can be computed (c.f. 
[H, m, m, [fii) from the evolution of A (Eq. (|B1|i 1. 
To do so, one needs to know how a stochastic differen- 
tial equation is written under a non-linear transformation 
since R and Q are nonlinear functions of the components 
of A. 

In general terms, let us call such a time dependent non- 
linear transformation ^{t, A) : A 1-^ ^{t, A) with compo- 
nents ^k, k g {1,2, ...,N}. Starting from the SDE of A 
(Eq. jE])), general formula [H, (52, [H, HI] give the new 
stochastic differential equations that governes ^, namely 

dt dA,j 2 dA.jdA 
" - - (B3) 



Dij'PgDpq'pg 



dt 



dA, 



dV 
'dt 



where the drift coefficients Ni are given by 
and the diffusion elements by 



3i? HfjAji — HfjAji 

3Q ^ H^jAjqAqi — HfjAjqAqi 



(B6) 



(B7) 



Afn = 2Tr [AA ' j + Q , 
M12 = A/21 = 2Tr f A^A^ 



1. 



A/22 = 2Tr [{A ' YA'j - 2Q' - -Tr (A^) 

Finaly, using again the general transformation Eq. (|B3[) . 
the probability current W of the joint probability 
V{Q*, R*) of the non-dimensional invariants Q* and R*, 
entering in the non-dimensional Fokker-Planck equation 



dt* {4, 



.w = o , 



(B8) 



is given by W = W(j]^.jf|^ + W^jjfj with (we recall that 
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drift 



w 



Q\R*\v{Q\R*)-{-^^ 



Q*,R* )V{Q*,R* 



(B9) 



W 



diflF 



We see from Eq. (jB9p that one needs to add another 
term W^jg- to the probabihty current W when deahng 
with a SDE. This term does not exist when deahng with 
a deterministic equation (Eq. ahhough, the gradi- 

ent of the forcing entering in the Navier-Stokes equations 
for A (Eq. ([T])) has been neglected. The total probabil- 
ity current W (Eq. (jB9[) ) was displayed in Fig. [Hh). 
We would like now to display separately the probability 
current coming from the drift terms VV^jj-jf^ and the one 
generated by the diffusion coefficients VV^jjff. We rep- 
resent in Fig. [To] the vector and streamline plots of the 
probability current for the model associated (a) to the 
drift coefficients W^jj-jf^ and (b) to the diffusion coeffi- 



cients Wjjff (Eq. (|B9p ). We recall that the total proba- 
bility current is displayed in Fig. [5]Jh). We see first that 
^drift ^^^^ ^diff ^'^'^ same order of magnitude. 

The current VV(^],j£^ associated to the deterministic part 
of the stochastic evolution (Eq. (jB4[) ) goes towards the 
origin in a rotating motion: the dynamics is decaying. 
To compensate for this decay, the current Wjjjff associ- 
ated to the stochastic forcing part of the evolution (Eq. 
(|B4[) ') points outwards away from the origin. The sum 
of these two, the total current displayed in Fig[5Uh), has 
a circular motion around the origin, consistent with a 
stationary process (i.e. dV{Q, R)/dt — — V.W ~ 0). 
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